Edit on GitHub

Examples

Reading Filterbanks

First, let's grab some filterbank data, we can find one from Berkeley's SETI Research Center. Julia's download function will grab the file from a URL and give us the filename back.

filename = download("https://github.com/UCBerkeleySETI/breakthrough/blob/master/GBT/filterbank_tutorial/blc04_guppi_57563_69862_HIP35136_0011.gpuspec.0002.fil?raw=true")
"/tmp/jl_5EO3ZV"

Now, we can read it into a Filterbank. By default, it will read all the time samples, we can however supply start and stop arguments if we want to look at a limited chunk. This is usefull as these files can get quite large.

fb = Filterbank(filename)
A Filterbank file with 272 time samples
----------
az_start            : 0.000000
foff                : -0.002861
nifs                : 1
za_start            : 0.000000
nbeams              : 1
machine_id          : 20
tsamp               : 1.073742
nbits               : 32
src_raj             : 71550.064000
src_dej             : 471420.040000
fch1                : 1501.463413
rawdatafile         : guppi_57563_69862_HIP35136_0011.0000.raw
ibeam               : 1
tstart              : 57563.808588
data_type           : 1
telescope_id        : 6
source_name         : HIP35136
nchans              : 65536

The data is stored as a DimArray from the DimensionalData package in the data attribute. This provides zero-cost abstractions for named indexing and fast index lookups. Our filterbank data has axes Freq and Ti, which we can use for indexing.

fb.data
272×65536 DimArray{Float32,2} with dimensions: 
  Ti Sampled 0.0:1.0737418239999998:290.984034304 ForwardOrdered Regular Points,
  Freq Sampled 1501.4634132385254:-0.00286102294921875:1313.9662742614746 ReverseOrdered Regular Points
 3.04948f9  3.14116f9  3.12580f9  …  8.60524f7  8.72117f7  8.30609f7
 3.00483f9  2.98634f9  2.97949f9     8.54793f7  8.44453f7  8.84643f7
 3.03344f9  3.00281f9  3.08042f9     8.42302f7  8.24331f7  8.63357f7
 3.01833f9  3.02201f9  3.03027f9     8.60705f7  8.34739f7  8.34594f7
 2.94578f9  2.99773f9  3.02747f9     8.04174f7  8.11727f7  8.16072f7
 2.92023f9  2.92796f9  2.91467f9  …  7.80949f7  7.99472f7  7.81455f7
 2.88401f9  2.95869f9  2.97999f9     7.91358f7  8.15653f7  7.97255f7
 ⋮                                ⋱                        ⋮
 2.95995f9  2.90506f9  3.01727f9  …  8.05276f7  8.04053f7  8.1892f7
 2.91456f9  2.92593f9  2.9851f9      7.91923f7  7.87905f7  7.93247f7
 2.89651f9  2.81591f9  2.85835f9     7.95392f7  8.06558f7  7.94178f7
 2.88585f9  2.86648f9  2.87025f9     8.20846f7  8.16997f7  8.20432f7
 2.90665f9  2.81179f9  2.91064f9     8.04896f7  7.90742f7  8.22847f7
 2.9013f9   2.93581f9  2.86071f9  …  8.0365f7   8.0742f7   8.21654f7
 2.94684f9  2.95507f9  2.94315f9     7.8024f7   7.79604f7  7.95808f7

Additionally, DimArrays have plot recipies that make visualization super easy. We can call Plots.jl's heatmap function to get a waterfall.

using Plots
fb.data |> heatmap

We can use this indexing to perform all sorts of analysis. Here we are using .. from IntervalSets.jl to create the range object.

fb.data[Freq = 1350..1450] |> heatmap

We can also index by time to get the nth integration

fb.data[Ti = 99] |> plot

We can also combine this indexing with Julia's builtin array operations such as looking at the time-averaged slice of spectrum. In this window, we get a nice look at the 21-centimeter line

using Statistics
dropdims(mean(fb.data,dims=Ti)[Freq = 1419..1422],dims=Ti) |> plot

Dedispersing

Let's look at some "realistic" data!

First we'll generate a fake dispersed pulse with a dispersion measure (DM) of 800 pc/cc over a frequency range of 1200-1500 MHz. We'll leave the other settings as is.

frb = pulse(800,1200,1500)
A Filterbank file with 2048 time samples
----------
tsamp               : 0.001000

Looking at the waterfall, we can see the frequency-integrated flux is spread out.

waterfall(frb)

We can apply the dedispersion transformation to get the real pulse shape

waterfall(dedisperse(frb,800))